module derivs
	use basicmodule
	implicit none
	
	contains
	
	function getdz(m,sig,xsi,y) result(val)		!derivative of 1+xsi*(y/sig-m)
		real	:: m,sig,xsi,y,val(3)
		val=[-xsi,-xsi*y/sig**2,y/sig-m]
	end function

	function getdw(m,sig,xsi,y) result(val) !derivative of (1+xsi*(y/sig-m))**(-1/xsi)
		real	:: m,sig,xsi,y,val(3)
		real	:: z,zx
		z=1+xsi*(y/sig-m)
		zx=z**(-1/xsi-1)
		val=[zx,zx*y/(sig**2),zx*(m*sig*xsi - xsi*y + sig*z*Log(z))/(sig*xsi**2)]
	end function
	
	function getd_logdens_base(msx,Y) result(val)
        real    :: Y(k),msx(3),m,sig,xsi,val(3),rk,zk, ri(k),zi(k),drk,dri(k)
        m=msx(1); sig=msx(2); xsi=msx(3)
		rk=y(k)/sig-m
		zk=1 + rk*xsi
		if(zk<0 .or. 1+xsi*(y(1)/sig-m)<0) then
			val=0
			return
		endif
		drk=(zk)**(-(1 + xsi)/xsi)
		ri=y/sig-m
		zi=1+ri*xsi
		dri=-((1 + xsi)/zi)
		val(1)=-drk-sum(dri)
		val(2)=(-drk*y(k)-sum(dri*y))/sig**2-k/sig
		val(3)=(-sum((ri*xsi*(1 + xsi))/zi) + sum(Log(zi))  + (rk*xsi - zk*Log(zk))/(zk)**((1 + xsi)/xsi))/xsi**2
	end function
	
	function getdlog_meanadj(Yk,msx) result(val)
        real    :: Yk,msx(3),m,sig,xsi,zk,rk,drk,val(3)
        m=msx(1); sig=msx(2); xsi=msx(3)
		rk=yk/sig-m
        zk=1+xsi*rk
        if(zk<0) then
           val=0
		else
			drk=((m + rk)*(-1 + xsi))/((1 + m + rk - m*xsi)*zk)
            val(1)=-drk+(1 - xsi)/(1 + m + rk - m*xsi)
			val(2)=-drk*yk/sig**2+1/sig
			val(3)=-((xsi*(-(rk*(1 + m + rk)) + xsi + rk*(2 + 2*m + rk)*xsi + rk*(1 - m + rk)*xsi**2) - (-1 + xsi)*(1 + m + rk - m*xsi)*zk*Log(zk))/((-1 + xsi)*xsi**2*(1 + m + rk - m*xsi)*zk))
		endif
	end function
	
	
	function getdlog_dens0(msx,Y) result(val)
        real    :: y(0:k),msx(3),val(3),ma
		val=getd_logdens_base(msx,Y(1:))
		ma=getmeanadj(msx,Y(k))
		val=val-2*(y(0)+ma)*getdlog_meanadj(Y(k),msx)*ma 
	end function
	
end module